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SUMMARY 


The central objective of this project has been the development of geostatistical methods for mapping elevation 
and ice surface characteristics from satellite radar altimeter (RA) and Synthetic Aperture Radar (SAR) data 
The main results are an Atlas of elevation maps of Antarctica, from GEOSAT RA data, and an Atlas from 
ERS-1 RA data, including a to ,al of about 200 maps with 3 km grid resolution. 

Maps and digital terrain models are applied to monitor and study changes in Antarctic ice streams and 
glaciers, including Lambert Glacier / Amery Ice Shelf, Mertz and Ninnis Glaciers, Jutulstraumen Glacier, 
Fimbul Ice Shelf, Slessor Glacier. Williamson Glacier and others. 
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(1) INTRODUCTION: MAPPING AND MODELING OF THE ANTARCTIC ICE SHEET 

The cryosphere plays a key role in the unstable global climate system. The polar ice caps and the Greenlandic 
inland ice shield are sensitive t.c changes in temperature (Huybrechts 1993). A collapse of the West Antarctic 
Ice Sheet could cause as much as 6 m sea-level rise (Bindschadler 1991). Discussion of instabilities of the 
Antarctic Ice Sheet was put forward as early as 1962 by J. Hollin. The mechanisms that may lead to ice-sheet 
collapse have been investigated and modeled in many studies, but are still a matter of debate (e.g.. Alley and 
Whillans 1991). Earlier work culminated in conclusions of catastrophic consequences (Hughes 1973; Mercer 
1978; Thomas 1977), whereas later modeling work showed that such catastrophic behavior is unlikely (Van 
der Veen 1985; Muszynski and Birchfield 1987). Scenarios for dynamic instabilities (ice-creep instabilities or 
even surging of the East Antarctic Ice Sheet) have been discussed (Clarke et al. 1977; Schubert and \uen 
1982). Huybrechts (1993) cone udes from a simulation that most of the variability in Antarctic ice mass and 
hence in sea level results from changes in the West Antarctic Ice Sheet, whereas the East Antarctic Ice Sheet 
seems to be robust to temperature changes. Changes in the East Antarctic Ice Sheet are also discussed in 
Colhoun (1991). The hypothesis that the stability of an ice sheet grounded below' sea level depends on the 
stability of its marine ice shelves (e.g., Mercer 1978; Thomas and Bentley 19/8, Lingle 1984) indicates a 
need to study the Antarctic ice shelves. That rapid retreat does occur in present times is documented by 
the examples of catastrophic ietreat. of Columbia Glacier, Alaska (Meier and Post 198/) and of bieak-up 
of Wordie Ice Shelf, Antarctic Peninsula (Vaughan 1993) (both located in warmer climates). Ice streams 
moving 10 to 100 times as fast as the adjacent, ice result in instability points in the dynamic system of an ice 
sheet. A prediction based on any model, however, can only be as good as the information on which the model 
is based. Many studies suffer from the fact that they are simulations lacking adequate data support. Satellite 
observations provide an efficient source of information for remote areas, and for large parts of Antarctica 
they are the best information presently available - once we understand how to use it right. One problem 
with investigations of the Antarctic ice mass is the lack of accurate topographic maps for large parts of the 
continent. 

The widely used Antarctic glaciological and geophysical folio edited by Drew-ry (1983) contains maps of a 
small scale only. A topographic map of the Filchner-Ronne-Schelfeis based on satellite images and ground- 
based geodetic surveys was recently published by Sievers et al. (1993). 

Most satellite payload yields images. Analysis of image data has many useful applications. Images of Antarc- 
tic ice streams have been compiled using AVHRR (Absolute Very High Resolution Radiometer) data from 
a NOAA (National Oceanic and Atmospheric Administration, USA) satellite (Bindschadler and Scambos 
1991) with a 1-km spatial resolution. Images of higher resolution are obtained by the Synthetic Aperture 
Radar (SAR) data, which have become available to the scientific community through ERS-1/2, JERS-1 and 
RADARSAT (ESA 1992a, b, 1993; Canadian Space Agency et al. 1994). A major difficulty with the analysis 
of SAR data is that quantitative analysis is not directly possible. One promising avenue in that direction 
is the application of interferometry, a technique that exploits the phase differences of two images, but at 
the same location, possible in the rare situation of very close repeat of the ground tracks (Goldstein et al. 
1993) and good correlation of the images to be compared (Zebker and Villasenor 1992). The best-known 
application is the extraction of the velocity of the ice (Goldstein et al. 1993). If no movement occurred 
and the environment did not change between the times of collection of the two images, it is possible to 
compute topography from pai -s of SAR images using interferometry. There is ongoing w'ork on construction 
of elevation maps from SAR stereo images, but that has yet to be completed. Examples of applications of 
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interferometry are restricted t.c date to the study of smaller regions, and SAR images can only be collected 
for 10 minutes per revolution. The technique is not suitable for mapping large areas of t he Antarctic ice. 
leaving ample necessity for altimetrv-based mapping. 

The best data source for topographic mapping from satellite is altimetry. The first satellite carrying an 
altimeter became operational in 1978 (SEAS AT). Together with data from the GEOSAT Geodetic Mission 
(1985-86) and the Exact Repeat Mission (1987-89) and data from ERS-1 (1992-96) and ERS-2. almost a 
20-year record of altimeter da? a is available. This makes altimeter data the type of data most suited for 
the study of changes on a regivinal or continental scale for length of record. One disadvantage of studying 
Antarctica by satellite data is t lat the orbital coverage of the previously mentioned satellites does not extend 
to the poles. 

Geostatistical analysis of satellite radar altimeter data may be utilized to construct maps of 3-km-by-3- 
km resolution of areas several 100 km large, which have a high accuracy (Herzfeld et al. 1993, 1994) 
Bamber (1994) produced a ma > of Antarctica (north of 82° S) from ERS-1 altimeter data with 20-km grids. 
Limitations of this map are tie lower resolution and the fact that the map is only reliable in areas with 
a slope of less than 0.65° (Bamber 1994). By total area most of Antarctica is flatter than 0.65° , but the 
steeper regions include the dynamically important ice streams and outlet glaciers. 

The geostatistical method (cf. Herzfeld et al. 1993) facilitates calculation of maps of higher accuracy and 
including steeper areas, but is computationally more intensive. The need for higher resolution is not well 
met if all of Antarctica is shown on one map sheet. An alternative is to construct an atlas, which in turn 
requires specific cartographic considerations. 

The central task of this projec . has been the calculation of an atlas of Antarctica, consisting of DTMs and 
maps with 3 km resolution, from GEOSAT and ERS-1 radar altimeter data, along with development of the 
necessary processing tools and geostatistical methods; and resulting in glaciologic applications. 


(2) REPORT: PREVIOUS WORK AND OBJECTIVES 

Work under this project built cn the development of geostatistical estimation (interpolation / extrapolation) 
methods and numerical implementation, specifically for the analysis of satellite radar altimeter data, resultant 
from my work under a previous NASA grant, for Lambert Glacier / Amery Ice Shelf (Herzfeld. Lingle and 
Lee, 1993, 1994). This metho l facilitates construction of digital terrain models with 3 km grid distance, 
high accuracy (50 m elevation accuracy on Lambert Glacier), and detection of the grounding line. Our result, 
of a 10 km advance of Lambert Glacier (by change of grounding line position) settled the at the time open 
question of advance or retreat of Lambert Glacier (based on several cross-over analyses). Our result could 
also be confirmed from cross-over analysis (Lingle et al., 1994). 

The original proposal NAGW-3790 (this project) had two central themes: (a) application of the geostatistical 
method to selected Antarctic glaciers and ice streams (Lambert Glacier, West Antarctic Ice Streams) and 
monitoring these glaciers; (b) development of a method to analyze SAR data. Following a request by NASA 
Polar Program manager Dr. Robert Thomas, the first objective was extended to evaluate all available radar 
altimeter data with the geostatistical method and produce maps for all of Antarctica, (north of the limit 
of satellite RA coverage, 72 . V S for GEOSAT, 82.1° S for ERS-1), and the second objective was largely 
dropped (but see section (5.1) on geostatistical classification for results). 
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Mapping ail of Antarctica with individual maps was a monumental task. 


(3.) GEOMATHEMATICAL TOOLS FOR RADAR ALTIMETER DATA ANALYSIS 

Geomathematical tools that were needed and developed for this project and are now available for other 
satellite data evaluation include: 

- geostatistical interpolation (with search routines adapted for RA tracks) 

- track-error correction routines 

- atlas-mapping scheme 

- TRANSVIEW tool 


(3.1) Data acquisition and correction 

A direct data connection between the Ice Sheet Altimetry Group at NASA GSFC (Dr. H.J. Zwally, Dr. J. 
DiMarzio and coworkers) and my group was set up for transfer of radar altimeter data. When ERS-1 data 
became available, our group was also instrumental in testing the necessary new correction algorithms and 
writing routines to identify bad tracks. Data processing by the Ice Sheet Altimetry Group includes: using 
the method of Martin and otlers (1983) for retracking, Goddard Earth Model (GEM) T2 orbits (Marsh 
and others 1983) for data reduction, and applying corrections for atmospheric effects and solid earth tides 
as described by Zwally and others (1983), slope corrections as described by Brenner and others (1983) and 
water-vapor corrections. After obtaining Ice Data Records (IDR) data sets, those points with retracked and 
slope-corrected data were retained. For each map sheet, a track plot is constructed to investigate coverage 
and ensure that coverage by retracked and slope-corrected data is sufficient. This was the case for all map 
sheets. After this processing, 11 bad” tracks with elevation (a) much lower than the surrounding area, or (b) 
of about constant small (50 rr ) difference to the surrounding area remained in several ERS-1 maps. An 
algorithm was developed to identify and remove these bad-track data. Elevation is given with reference to 
the WGS84 ellipsoid. 

(3.2) Atlas mapping for continent-wide coverage 
(3.2.1) Atlas concept 

Rather than inverting all the Antarctic radar altimeter data onto a grid to produce a single map covering 
Antarctica (with, of course, a hole for the area south of the limit of radar altimetry coverage), we use an 
atlas mapping scheme. This improves resolution, facilitates mapping of detailed structures, and reduces 
distortion due to cartographic projection, which is particularly severe for high latitudes and large areas in 
most algorithms. 

An atlas in the sense of diffeiential analysis (Holmann and Rummler, 1972, p.63) is a set of maps that 
(i) covers a given area completely (that is, each point in the area is contained in at least one map); («') 
projections restricted to areas that appear on two (adjacent) maps (subsets of two maps) are identical on the 
intersection. In cartography, only property (/j is required for an atlas, and the neighbourhood relationships 
need to be matched between sheets. 



A useful projection for mapping at high latitude and in atlas form is the Universal Transverse Mercator 
Projection (UTM) (Snyder, 1987), which results in an orthogonal coordinate system with coordinates in 
meters. Sufficient overlap of acijacent sheets is convenient for the user of the atlas and necessary to ensure 
that each point of Antarctica s contained in at least one map despite of the distortion of the map edges 
introduced by the projection algorithm. 


(3.2.2) TRANSVIEW tool 

One of the oldest problems in mapping the Earth is the definition of projections of the Earth’s surface onto 
a two-dimensional map sheet. For mapping purposes, the Geoid is commonly approximated by a sphere 
or ellipsoid. Desirable properties of map projections are conservation of area (equal-area projection), of 
distances (equal-distance projection), of angles (equal-angle or conformal projection), which are mutually 
exclusive when mapping oil a plane, and projection to a rectangular coordinate system (for examples see 
Hake, 1982; Snyder. 1987). Because it is not possible to satisfy all of these conditions, some projections have 
been defined that do not fulfill any conditions exactly, but a combination of them approximately (Hake, 1982: 
Snyder, 1987). For series of topographic maps in countries with a long tradition in mapping, algorithms 
have been designed to construct maps constituting an atlas. 

For mapping Antarctica, however, we had to design a much-needed tool to convert, match and visualize 
UTM coordinates and geographic coordinates, to satisfy the Atlas mapping conditions. 

From the viewpoint of interpo ation of irregularly distributed data onto a regular grid, an orthogonal co- 
ordinate system facilitates the algorithm and saves computation time. The latter is especially important if 
distance-dependent measures are used, such as in inverse-distance weighting or in geostatistical methods. 

Distance may be calculated on ' he sphere or on the ellipsoid (cf. Moritz, 1980; Torge, 1980), but this requires 
transformations at each step cf the interpolation algorithm which usually is dependent on the number of 
points squared. In comparison : the number of essential operations for coordinate transformation depends 
only linearly on the number of points. Methods involving the covariance function or the variogram (kriging. 
least-squares prediction; cf. F erzfeld, 1992) would require estimation of the structure function over the 
ellipsoid which would be troublesome. It is apparent that coordinate systems with orthogonal coordinates 
in meter units are thus extremely convenient for interpolation purposes. 

Common practice is not to change coordinate systems, but to simply use geographic coordinates, which is 
unproblematic for small areas. For large areas, neglection of the coordinate transformation results in a severe 
distortion of the spatial structure in the data. (Recall that 1 0 latitude is always about 111 km, but 1 0 
longitude is cosine of latitude times 111 km; so, at 60 0 North/South it is only 0.5 times 1 1 1 km or 55.5 km.) 

The distortion is particularly severe for mapping at high latitude. The Arctic and Antarctic are usually 
treated separately in one map using the polar stereographic projection. Typically, such maps of polar regions 
are at a small scale and do not show much detail. The importance of the polar system in the Earth’s global 
systems and its role in ’global change’ have become increasingly recognized. A useful projection algorithm 
for mapping at all latitudes is the Transverse Mercator projection. A Mercator projection is defined by a 
cylinder that is tangent to the Earth and a mapping to orthogonal coordinates. For the (common) Mercator 
projection, the tangent circle is the Equator, for the Transverse Mercator projection, the tangent circle 
is a meridian (called the cential meridian of the projection). The advantage of the Transverse Mercator 



projection is that ali latitudes are mapped with the same distortion. The disadvantage is that areas far 
away from the central meridian are strangely distorted. The solution provided by the UTM system is to 
rotate the cylinder around the Earth in steps of 6 0 . Zone 1 corresponds to central meridian 177 0 W. The 
standard central meridians arc at 3 0 (for 0 ° - 6 0 ), 9 0 (for 6 0 - 12 0 ), 15 0 (for 12 y - 18 0 ) etc. The 
central meridian is projected with a factor of 0.9996. lines of true scale are approximately parallel and lie 
approximately 180 km east and west of the central meridian. The border meridians are projected slightly 
lengthened, for example at 50 ° latitude with a factor of 1.00015. The projection is defined everywhere 
except 90 0 away from the central meridian. In the UTM scheme, the projection is chosen such that the 
central meridian is mapped tc East coordinate 500,000, units are in meters; the North coordinate along 
the central meridian is in meters from the equator (along the ellipsoid). UTM is conformal and close to an 
equal-area projection, it has orthogonal coordinates in meters. The UTM projection, for instance, does not 
satisfy condition (ii) of an atlas, if two adjacent maps belong to two different central meridians: however, 
this defect may be compensated for by overlaps between neighboring map sheets. 

The objective of our program TRANSVIEW is to provide a tool to calculate and visualize 

(a) the shape of a map that, is rectangular in geographic coordinates when transformed to UTM coordinates 

(b) the amount of distortion for any map sheet on the Earth 

(c) the largest map that is rectangular in UTM coordinates and inscribed in a given map that is rectangular 
in geographic coordinates 

(c) the overlap necessary to map a large area in individual sheets using the UTM projection, and 

(d) to provide a system that works also for Antarctica. 

TRANS VIEW works for any r 'ctangular area on the Earth. The only restriction is that the area needs to 
be located entirely on the Northern or on the Southern hemisphere. The UTM projection is defined relative 
to an appropriate central meridian (3° W, 9° W, 15° W, etc., uneven multiples of 3° ). For map areas that 
do not contain a central meridian of the UTM projection, an appropriate meridian needs to be determined 
by the user of the program. The latter is of particular importance for mapping small areas. The program 
discerns the location of the map relative to the central meridian and automatically selects the appropriate 
case for the transformation (se * Figure 1). All possible cases are given in Figure 1. The algorithm is given 
in Herzfeld et al. (1999) and available from the website of the International Association for Mathematical 
Geology (http://www.iamg.org/). An example of transformation and visualization is given in Figure 2. 


(3.2.3) Result: Definition c f map sheets for the Antarctic Atlas Projects 

GEOS AT ATLAS: For the GEOSAT Atlas, the rows and sheet tiling are defined as follows: rows: 72.1 0 - 
67 0 ; 68 0 - 63 0 ; 64 0 - 60 0 . Maps in row 72.1 0 - 67 0 are 546 km (183 gridnodes) E-W and 543 km (182 
gridnodes) N-S. Maps in row 6S 0 - 63 0 are 666 km (223 gridnodes) E-W and 531 km (178 gridnodes) N-S. 
Latitude 72.1 0 South marks the poleward limit of coverage of altimeter data from the Seasat and Geosat 
satellites. The size of sheets for the Antarctic atlas is defined as follows: 16 0 longitude span, 2 0 overlap on 
each side, 12 ° offset from one map to the next, and 1 0 overlap at top and bottom. 

ERS-1 ATLAS: Map sheets of the ERS atlases are designed to match those of the GEOSAT atlas. The 
Antarctic is divided into rows of map sheets (63° - 68° S, 67° - 72.1° S, 71° - 77° S, 75° - 80° S, 78° - 81.5° 
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cose 2, southern hemisphere, m69e57 73n72l -67(a) cose 2. northern hemisphere, m69#57 73n721-67(n) 




cose 3, southern hemisphere, m69e61 77n721-67(s) case 3. northern hemisphere. m69e61 77n72l-67(n) 




case 4, southern hemisphere, m69e65 81n721-67(s) cose 4, northern hemisphere, m69e65 B1n721-67(n) 




case 5, southern hemisphere. m69e73 89n721-67(s) case 5, northern hemisphere, m69e73 89n721-67(n) 




Fig. 1. Cases of possible location of ar area that is rectangular in geographic coordinates relative to central meridian. Program 
TRANSVIEW distinguishes automatically among these cases. Plots in UTM coordinates; solid circles: nominal map area that is 
rectangular in geographic coordinates; open circles: inscribed map area that is rectangular in UTM coordinates; vertical line: cen- 
tral meridian of the UTM projection {mapped to 500,000). 
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ANTARCTIC MAPPING PROJECT 


map view geo to utm coord trafo m69e61-77n721-67.dat (I6deg) 



Fig* 2, Example of transformation and visualization using TRANSVIEW, for area 61°-77° E/67°-72.1° S, using central meridian 
69° E (UTM zone 42) and map size of 16° east-west and 5.1° north-south. Curvilinear rectangle marked by filled circles delineates 
map edge of area originally rectangular in geographic coordinates after transformation to UTM system. Straight rectangle marked 
by empty circles delineates edges of inscribed UTM map. Annotated coordinates in brackets are backtransformed geographic coor- 
dinates; eastlap, toplap, botlap are minimal overlaps with adjacent sheets in degrees required for gap-free coverage in east-west, 
north and south direction, respectively 




S). In rows 63° - 68° S and 67°' - 72. 1^ S maps are of 16 degrees longitude nominal size and overlap 2 degrees 
on each side, so each map is cffset against the next one by 12 degrees longitude. In rows 7T J - 77° S. 75° 
- SO 0 ' S, and 78° - 81.5° S maps are of 36 degrees longitude nominal size, overlap is 6 degrees on each side, 
offset of two adjacent maps is 24 degrees longitude. The map names (e.g. m69e61-77n67-721) give central 
meridian (69° ) and extent of the nominal map area (61° to 77° E, 72.1° to 67° S (the minus in the north 
coordinate was dropped in th • file names for ease of reading)). The true map area is then the maximal 
area contained in the nominal map area with the property that it has straight edges in the UTM coordinate 
system and the grid coordinates are multiples of 3000 (Herzfeld and Matassa, 1999). All maps in one row- 
have the same size and grid coordinates. 

(3.3) Geostatistical methods for elevation mapping 

Interpolation of corrected altimeter elevation data is carried out using geostatistical methods, in this appli- 
cation, ordinary kriging. ’’Kriging” comprises a family of interpolation/extrapolation methods based on the 
least-squares optimization principle and usually introduced in a probabilistic framework (Matheron, 1963: 
Journel and Huijbregts, 1978) However, it can be considered as a numerical interpolation technique only 
(Herzfeld, 1992a). Kriging consists of two steps: 1. analysis of spatial structures (variography), and 2. 
estimation, interpolation and extrapolation (kriging). 

(3.3.1) Effect of variogram models on elevation mapping 

In the first step, a variogram is calculated according to 

7 (h) = A £ [*(*,) - *(*,- + h)f~ (1) 

i = 1 

where r(x t -), z(xi + h) are measurements at locations Xi , x, -f h, respectively, inside a region D , and n is the 
number of pairs separated by the vector h. The residual variogram is 

res{h) = 7 (h) - ^m(h ) 2 ( 2 ) 

where 

1 n 

™(h) = - X! [-(*•) - -(• r < + *)] (3) 

U 1 = 1 

is the drift component. The experimental variograms are calculated in distance classes from the measure- 
ments, in our case from the ra lar altimeter (RA) data, then are fitted by analytical variogram models. A 
variogram model describes the type of transition from the strong covariation between closely neigbouring 
samples to weaker covariation of samples farther apart. The variogram model is characterized by its func- 
tion type, which has to meet he positive-definiteness criterion to ensure existence and uniqueness of the 
solution of the kriging system. For mapping altimeter data of ice surfaces, a Gaussian variogram model with 
a relatively high nugget effect is used, because 

(1) the ice surface at 3 km resolution is smooth, and the Gaussian model is the model of a mean-square 
differential process, 

(2) altimeter data at this resolution have a high signal-to-noise ratio, and the spacing of tracks influences 
the variogram. 
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The function of the Gaussian variogram model is 


7 (/i) — C'u + Ch (1 — exp( -•)) (4). 

a~ 


The nugget effect, Co, is the residual variance of resampling in the same location for altimeter data, it 
also contains contributions caused by surface morphology at a resolution higher than the resolution of the 
altimeter footprint (called the support). The total sill, Co + Ci. is close to the total variance of the data. The 
concept of a regionalized variable, which is fundamental in geostatistics, perceives that closely neighbouring 
measurements have a higher covariance than measurements spaced farther apart. For distances increasing 
from the support to the range parameter a in equation (4), the variogram increases and the regionalization 
effect decreases; for distances beyond the range, the regionalized variable theoretically behaves like a random 
variable (probabilistically speaking), and data with a distance larger or equal to the range all have the same 
influence on the interpolated value (numerically speaking). 

For GEOS AT data 

C’o = 250m 2 

C\ = 343m 2 (5) 

a = 18000m; 


for ERS-1 data 

C’o — 43m 2 

Ci = 18 m 2 (6) 

a — 16000m. 


The noisiness of the data is captured by the ratio no of the nugget effect Co to the total sill Co T C\ in the 
following equation: 


no = 


Co 

Co T Ci 


( 7 ). 


The nugget effect ratio no is 0.422 for GEOSAT data and 0.705 for ERS-1 data, so ERS-1 data are noisier. 
The variograms were fitted to data from an ice stream around the grounding zone (Lambert Glacier). To 
avoid edge effects between maps of the same atlas, the same variogram is applied throughout any one atlas. 

The effect of using a variograir with a higher nugget effect ratio is that contour lines are smoother. This is 
correct, because noisiness in the data does not warrant to pick up too many details everywhere. 

Selection of the variogram, however, does not only depend on the data analysis, but also on the mapping 
purpose. For monitoring of cfanges from one year to another, the same variogram needs to be used for 
both years, to facilitate companson. This has been applied in the study of Lambert Glacier/Amery Ice Shelf 
(Herzfeld and others, 1997). 

To improve the representation of ice surface structures in terrain models of individual glaciers, regional 
variations of the variogram need to be taken into account. To this extent, classes of ice with homogeneous 
surface properties such as floating ice tongues, grounded glacier ice, steep margins, mountainous terrain, 
and almost flat inland ice are characterized by specific variogram functions, which may then be used in the 
interpolation. 
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(3.3.2) Kriging 


The kriging method called ’‘ordinary kriging” (universal kriging of order 0) is better suited for interpolation 
of RA data than universal knging of a higher order, because the drift parts modelled by a polynomial 
component in the higher order universal kriging methods is likely to create artefacts in the gaps. 

In ordinary kriging, the value ::y = c(j?o) at a node £o is estimated by 

71 

Zq = QtjZj with a, e IR (8) 

*=i 

with data — £(. l * i ) at locations Xi [i — 1, n) in a neighbourhood of the grid node ^o- In the Antarctic atlas 
mapping, the 4 nearest neighbours in each quadrant were used in the interpolation of any given grid node. 
The coefficients are determined such that the estimation error has minimum variance and the estimation is 
unbiased, which requires a condition 

n 

= 1 (9)- 

1 = 1 

A solution of the kriging system is obtained using the variograrn model specified earlier. (For a derivation of 
the kriging system, see Herzfeld 1992a,b.) This is carried out for every grid node in the map area. A 3 km 
grid size is chosen so that the resultant maps or terrain models can be used in glaciodynamic investigations. 


(4.) GLACIOLOGIC RESULTS AND APPLICATIONS 
(4.1) GEOSAT and ERS-1 Atlas of Antarctica 

The GEOSAT Atlas was produced from 1985-86 GM data, the ERS-1 Atlas from 1995 GM data. Work on 
the ERS-1 Atlas was continued after completion of the work under this grant. Both atlases contain a total 
of over 200 individual maps (some are given in Figures 3-8) (see also Herzfeld and Matassa, 1999; Herzfeld 
et al., 2000a). 

The Atlas DTMs are available through the website of the National Snow and Ice Data Center 
(http:/ /www. nsidc.org/). 

Quality of maps: Maps are available as either contoured maps or digital terrain models. For the study of 
individual glaciers, it should noted that much more information is contained in the Atlas DTMs than is 
visible in the Atlas maps. For example, Williamson Glacier: Williamson Glacier flows into Moscov University 
Ice Shelf at a location east of Law Dome, as seen on the ERS-1 atlas map ml 17e 109-1 25n63-68 ‘‘Sabrina 
Coast” (Fig. 5a in Herzfeld el al., 2000a). The details in the digital terrain model may be investigated, 
as shown in the map of Williamson Glacier in Figure 5b in Herzfeld et al. (2000a). The resolution of the 
detail map is the same (3 km 'rid) as that of the entire atlas. This demonstrates that more information is 
contained in the digital terrain model than is obvious from a plot of the maps at atlas scale. The models 
may be used for geophysical study, or for monitoring smaller outlet glaciers of the Antarctic Inland Ice. Not 
much information other than satellite data is available for this part of East Antarctica. 

So, there are numerous glaciol >gic applications possible. A few are given in the sequel. 


10 



north (km UTM) 










north (km UTM) 





north (km UTM) 






north (km UTM) 


Lambert Glacier - ERS1 DATA, 1995 


-7500 


-7600 


-7700 


“7800 


-7900 



300 400 500 600 700 

east (km UTM) 

e61 -77n67-721 , WGS84, Gaussian variog., central mer. 69, slope 
corrected, sccle 1:5000000, 970727 





5600 

5400 

5200 

5000 

4800 

4600 

4400 

4200 

4000 

3800 

3600 

3400 

3200 

3000 

2800 

2600 

2400 

2200 

2000 

1800 

1600 

1400 

1200 

1000 

800 

600 

400 

200 

140 

130 

120 

110 

100 

90 

80 

70 

60 

50 

40 






north (km UTM) 


-8500 


-8600 


-8700 


-8800 



Executive Committee Range - ERS1 DATA, 1995 



200 300 400 500 600 700 800 

east (km UTM) 


e21 9-255n75-80, WG584, Gaussian variog., central mer. 237, slope corrected, scale 
1:5000000, 971 105 


5600 

5400 

5200 

5000 

4800 

4600 

4400 

4200 

4000 

3800 

3600 

3400 

3200 

3000 

2800 

2600 

2400 

2200 

2000 

1800 

1600 

1400 

1200 

1000 

800 

600 

400 

200 

140 

130 

120 

110 

100 

90 

80 

70 

60 

50 

40 



(4.2) Presentation of ice surfaces witli reference to Geoid models 


In the previous sections, elevations referenced to the WGS84 ellipsoid have been employed. When studying 
changes, the reference surface is unimportant. Comparison with published studies and field work may be 
facilitated by maps referenced i.o a geoid, although for the Lambert Glacier area elevations are rarely given 
in the literature. Referencing to a geoid, however, involves geodetic problems of geoid determination, which 
is difficult in areas of poor availability of gravity data. Mathematical methods for geoid determination are 
described in Moritz (1984), Ts< herning (1984), and Engels and others (1993). 

We used two recently published geoids, ” Goddard Earth Model” GEM-T3 (Lerch and others, 1992) and 
’’Ohio State University Model’' 0SU91A (Rapp, 1992, 1994), kindly made available to us by NASA GSFC 
and by R. H. Rapp (pers. comm., Feb. 1995). Lambert Glacier maps referenced to GEM-T3 and 0SU91A 
are given in Figures 3F and 3G in Herzfeld et al. (1997), respectively. Comparison of the two geoid models in 
the area of our maps shows significant differences: The GEM-T3 surface gradually slopes in a northeasterly 
direction, the 0SU91A surface has a more complex topography with a hill (Herzfeld et al.. 1997, 1998). 
According to R. H. Rapp (pers comm., Feb. 1995), both models are based on the same set of data of only 
fair reliability (mean anomalies in 1° by \° cells), therefore differences between the two models should 
be attributed to philosophical differences in the interpolation method rather than warrant a geophysical 
interpretation. For the OSU91A model, global average in estimation of the total undulation error is 0.57 m. 
but 2 m for land area with no surface gravity data (Rapp, 1992). 

On the other hand, our ellipsoid-referenced topographic maps can be used for terrain correction in the inverse 
gravimetric problem, and thus may facilitate improvements of the geoid in this poorly constrained region. 
For purposes of field study ar d comparison to the literature both geoid-referenced maps are sufficiently 
accurate, as well as the WGS ■'T -referenced maps when a constant of about 20 m is subtracted from the 
elevation values (the approximate difference between ellipsoid and geoid in the map area). 


(4.3) Monitoring Lambert Glacier / Amery Ice Shelf system 

Lambert Glacier / Amery Ice Shelf is the Antarctic ice-stream/ice-shelf system with the longest record of 
satellite RA data coverage, starting with SEASAT (1978), simply because it is located largely north of 72.1° 
S. By analyzing time series of satellite-altimetry- based digital terrain models, we studied changes in Lambert 
Glacier. Data for austral winters 1978, 1987, 1988, 1989 were used and various types of correction studied 
(total of 17 maps) (Herzfeld ft al., 1997). When ERS-1 data became available, monitoring of Lambert 
Glacier was continued with 1992 maps and a 1995 map (Matassa et al., 1996; Herzfeld and Matassa, 1999). 

A detailed discussion of sensitivity of the results to slope correction, frequency filtering, and mathematical 
method and statistical algorithms for calculation of changes is given in Herzfeld et al. (1997). As a result 
it was found that data that are neither slope-corrected nor frequency-filtered are best suited for monitoring 
changes. 

The purpose of calculating elevation changes is to answer the following glaciological questions: 

(GQ1) Time distribution of changes: How is the general advance of the grounding line between Seasat 
and Geosat times (Herzfeld, L ngle, and Lee, 1994) and overall increase in surface elevation, deduced from 
crossover analyses ( Lingle and others, 1994) distributed in time? 
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A general increase in elevation for the entire glacier/ice-shelf system is derived from averaging elevations ot 
the DTMs. This is consistent, with results from crossover analysis (Lingle and others. 1994) and the advance 
of the grounding line deduced from DTMs (Herzfeld, Lingle. and Lee. 1994). The absolute differences 
between Seasat and Geosat ck La, however, are lower than the fluct uations between the three Geosat years 
1987, 1988, and 1989, with 1988 values either too high or too low to match the trend (see Table 1). These 
may be due to interannual variation, or more likely, to undetected problems with the altimeter, because 
apparent interannual elevation changes indicated by the 1988 figures are rather large. With available data, 
however, we have no means to exclude a sudden change. 

(GQ2) Location breakdown: Is the glacier homogeneously increasing in elevation ? 

In summary, a significant increase in elevation is observed for the Amery Ice Shelf and the grounding zone 
(for all outline locations and both averaging algorithms tested), whereas average surface elevations of the 
rougher and topographically more varied grounded Lambert Glacier may have decreased (but. the decrease 
is not sufficiently supported by the data). 

It is concluded that the glacier advanced about 10-12 km between Seasat and Geosat ERM times (taken 
for all or each of 1987, 1988, 1989), and that the apparent interannual variability between 1987. 1988, and 
1989 is high relative to the change observed for the 10-to- 12-year interval. 

Partington and others (1987) identify the position of several points along the grounding line. Their map. 
however, does not show any coordinates nor scale, which makes an exact comparison impossible. Their 
northernmost grounding point is in the same general area (south of Beaver Lake and Jetty Peninsula) as 
the grounding zone on our Seasat-data-based map (Fig. 3 A in Herzfeld et al.. 1997), the other points are 
up to about 80 km farther soui-h. Budd, Corry, and Jacka (1982) identify a grounding point from a break in 
surface slope at LTTM 485,600 F/-7902.900S (3 km south of their point T4), this grounding point is located 
in the southern part of the grounding zone identified on our Seasat-data-derived map (Fig. 3A in Herzfeld et 
al., 1997). This comparison might indicate a small advance bet ween the time of the Australian expeditions 
(1968-1971) and 1978. Brooks and others (1983) find a break-in-slope 43 km south of the grounding point 
of Budd, Corry, and Jacka (1982); however, Brooks and others (1983) remark that their altimeter data 
processing does not allow to identify the grounding line. According to Partington and others (1987) most 
estimates in the literature indicate a net accumulation, except McIntyre (1985). How much of the observed 
elevation changes and the advance of the grounding line is due to mass increase cannot be determined from 
the available data. A more detailed survey of the velocity field of Lambert Glacier/ Amery Ice Shelf may 
contribute to answering related questions. 

The dynamics of the Amery Ice Shelf and the ice streams that feed into it may be complex. The fast-flowing 
ice discharges at about 1.2 kn/’yr into the shelf, only about one third of the discharge is contributed by 
(lower) Lambert Glacier (Budd, Corry, and Jacka, 1982). Wellmann (1982) finds geomorphologic evidence 
that Fisher Glacier, one other tributary of the Amery Ice Shelf, may be surging. Budd (1966) reports a 
50 km retreat of the ice shelf front between 1955 and 1965. Consequently, trends in elevation changes need 
not necessarily coincide for th ■ lower Lambert Glacier and the Amery Ice Shelf. This same observation is 
made in our analysis. 

Brooks and others (1983) conclude from a comparison of Seasat data and data published by Budd, Corry. 
and Jacka (1982) that Lambei t Glacier may have retreated between 1968-1971 and 1978. while Charvbdis 
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Glacier (another tributary) may have surged. However, a visual comparison between the Seasat-altimetry- 
derived map by Brooks and o’ hers (1983) and our Seasat-based map (Fig. 3A in Herzfeld et al.. 1997). 
using the *20 m-mean offset, between the Geoid and the WGS84 ellipsoid (cf. Figs. 3E and 3F in Herzfeld 
et ah, 1997), shows that the shapes of the contour lines match well in the floating part and that absolute 
differences are within 5 meters, but differences in morphology exist in the grounded part of Lambert Glacier 
(except for a general high on ti e western side of Lambert Glacier, and a low on the eastern side, common to 
both evaluations). The differences are a consequence of the different approaches to evaluation. Brooks and 
others (1983) used an inverse-distance gridding algorithm and contoured by hand, based on elevations from 
orbits adjusted into a common ocean surface. These results indicate that some of the differences in different 
studies may be due to different processing algorithms and data references and do not permit a glaeiodynarnic 
interpretation. 

Reasons for apparent rapid changes are difficult to determine from remote sensing information only. There 
is a possibility that a change in the dynamic system of the glaciers feeding into Amery Ice Shelf occured. 
Thickening of the lower part of the system and thinning of the upper part, as indicated by the spatial 
breakdown in our analysis, is typical for surges. The hypothesis of a dynamical event is mentioned here 
for sake of completeness, as the evidence from altimeter data alone is not sufficient for such a conjecture. 
The hypothesis is found in the literature (Wellmann, 1982: Brooks and others. 1983), but largely without 
supporting data. 

(4.4) Investigations of other glaciers 

Other Antarctic glaciers and areas studied and described in some detail include: 

- Mertz and Ninnis Glacier To agues (Herzfeld and Matassa, 1997) 

- Riiser-Larsen Peninsula (Herzfeld and Matassa, 1999) 

- Prince Olav Coast (Herzfeld and Matassa, 1999) 

- Mawson Coast West (Herzfeld and Matassa, 1999) 

- Ingrid Christensen Coast (Herzfeld and Matassa, 1999) 

- Pennell Coast (Herzfeld and Matassa, 1999) 

- Napier Mountains (Herzfeld and Matassa, 1999) 

- Knox Coast (Herzfeld and Matassa, 1999) 

- Sabrina Coast (Herzfeld and Matassa, 1999) 

- Graham Land, Antarctic Peninsula (Herzfeld and Matassa, 1999) 

- Slessor Glacier (Herzfeld et al., 2000a) 

- Fimbul Ice Shelf, Jutulstraur ten Glacier (Herzfeld et al., 2000a) 

- Williamson Glacier (Herzfeld et ah, 2000a) 
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(5.) WORK ON SAR DATA AND GEOSTATISTIC AL SURFACE CLASSIFICATION 

(5.1) Geostatistical Surface Classification 

The difference between geosta istical interpolation and geostatistical classification is explained in Herzfeld 
(1999). As a consequence of the shift in focus of the proposal from classification method development to 
Antarctic-wide RA data analysis, and due to the delay in availability of RADARS AT SAR data, not much 
effort could be expended under this project to the development of the geostatistical surface classification 
method for SAR data. Nevertheless, I give a brief summary of results on geostatistical classification and 
include those publications, in which the partial support through this grant is acknowledged. Through some 
of my other work, the principles of geostatistical classification have by now advanced to a considerable level, 
and some roots of these ideas go back to work under NASA grant NAGW-3790/NAG5-61 14. 

Geostatistical estimation (interpolation/extrapolation) is a known method, adapted by myself to the evalu- 
ation of RA data. Geostatistical classification now summarizes a suite of methods developed by myself for 
various classification projects in geophysics and glaciology. 

While interpolation utilizes the primary information in the data, a geostatistical surface classification method 
is developed to derive secondary information from elevation and backscatter data. Based on quantitative 
properties of the variogram, elements of surface structures are used for mapping and segmentation of an area 
into provinces homogeneous in surface characteristics. 

A critical issue in the analysis of satellite Synthetic Aperture Radar (SAR) data is the availability of ground 
truthing to distinguish between intensity variations caused by subscale-resolution geophysical variability and 
noise, and to determine small-scale sources of variations in backscattering. During the 1993-1995 surge of 
Bering Glacier, Alaska, GPS-1 seated video data were collected from small aircraft and analyzed system- 
atically with the geostatistical ice surface classification system ICECLASS. The objectives are to (a) help 
understand the relationships between ice velocities, surface strain states, and progression of deformation 
processes during the surge, and (b) provide a technique for surface classification based on video data, SAR 
data, or image data in genera!. 

For principles of surface classification, see Herzfeld (1998) (also images of Bering Glacier during the surge). 
A connectionist-geostatistical system for classification is given in Herzfeld et al. (1996). Application of 
classification to SAR data is the objective of Herzfeld et al. (2000b). 

(5.2) Mertz and Ninnis G acier Tongues — Comparison of results from SAR and RA data 
analysis 

A comparison of quantitative information obtained from SAR and RA data gave surprising results for Mertz 
and Ninnis Glaciers, Antarctica (Herzfeld and Matassa, 1997). 

Mertz and Ninnis Glaciers are both glaciers with long tongues that extend into the ocean and fluctuate 
considerably in length. Results from the GEOSAT Antarctic Atlas DTMs were compared to results from 
SAR data (Wendler et ah, 1996). Mertz and Ninnis Glacier maps were expanded from Atlas maps, and 
details of slope and length of I he tongue were measured. The detail maps showed a surprising amount of 
detail (which was later confirm- d by geologists working in the area; G. Kleinschmidt, pers. communication). 

Mertz Glacier is at UTM 360,010 - 440,000 E / 7,525,000 - 7,440,000 S. Its drainage is a broad valley that is 
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distinguishable at the southerr map boundary. The sides of the western valley are transected by erosional 
features. At its head, the glacier is fed by a narrow valley, approximately 4 km wide. The distance from the 
80 m contour line at the head o the ice front is 85 km. In the eastern arm of Mertz Glacier there is a 20 m 
overdeepening, centered at 415.000 / -7510. This is located below the steepest part of the valley walls. This 
is not an interpolation error, b *cause the contours are really smooth (see Figures 9-11). 

Mertz Glacier is about 27 km wide. The glacier tongue appears to extend about 45 km seaward of the 
coastline (see discussion below). The grounding line of Mertz Glacier is probably located in the vicinity of 
the 40-50 m contour. The 60 m contour still exhibits the indentation of the valley that continues subglacially 
from the valley leading to the head of the glacier, while the 50 m contour does not. The 40 m and 50 m 
contours show the signs of the eastern side valley. At 30 m, the tongue is definitely floating. Slopes to 2° 
were calculated accurately. 

Ninnis Glacier is at UTM 500.000- 540,000 E / 7,550,000- 7,580,000 S. The glacier tongue appears to extend 
about 10-15 km seaward from the coastline (the coastline being identified by the steep gradient of contours). 
Ninnis Glacier lies below a steep cliff, at 90 m above WGS 84 and lower. The entire extension from the 
break in slope at the foot of the cliff to the 0 m contour line is 20 km. [Notice that the 0 m contour does 
not coincide with sea level, because the reference level is the ellipsoid.] Ninnis Glacier is about 35 km wide. 
Ninnis Glacier does not really have a ” tongue*' anymore (as opposed to 1913) as noted by Wendler et al. 
(1996). 


Discussion on accuracy and reliability. Mapping of atlas type and careful analysis of the spatial structure 
and the distribution of noise levels in radar altimeter data allows us to extend the limits of use of altimeter 
data for mapping. Bamber (1994) produced a map of Antarctica (from ERS- 1 altimeter data) with 20-km 
grids, reliable only in areas with a slope of less than 0.65° according to the author. The drainages of Mertz 
Glacier and the Mertz and Ninnis Glacier area are considerably steeper than that. Accuracy depends on 
topographic relief, as calculated in Herzfeld et al. (1993. 1994), with submeter accuracy on ice streams. 
However, the accuracy of a kriged map in areas of high topographic relief refers to the pointwise error rather 
than to the error of the ”mean ! surface elevation mapped. The "pointwise error” is the probable difference 
between a point on the map (estimated surface elevation) and a radar- altimetry-derived surface elevation: it 
depends on the averaging process of kriging. Shape and elevation of the surface at the Mertz Glacier drainage 
and on the glacier itself, however, are more accurate than inferred from the noise calculation, which is best 
conceived from inspection of the maps. The contour lines are smooth on the resolution of the 3-km grid. 
Smooth contour lines indicate a continuous or differentiable surface function with low error; little islands 
and edging contour lines indiec.te a rougher surface function with higher errors. That means the maps are 
reliable also in areas of steep terrain. An area with high relief is located on the western side of Mertz Glacier. 
Notice that the grid spacing of the subarea enlargements is the same 3 km as for the large map. and that a 
wealth of information becomes only available in the enlargement. Mertz Glacier Tongue appears to extend 
about 40 km seaward of the coastline, Ninnis Glacier Tongue about 20 km. 

Neither satellite altimetry nor vriging are tools designed to track the location of an ice cliff. Because of the 
effect described in Thomas and others (1983) and Partington and others (1987) (snagging of altimeter) the 
location of the ice edge is systematically wrong, differently so in descending and ascending orbits. Thomas 
and others (1983) attempted to trace the ice edge from altimeter data designing a technique not used here. 
In the Ice Data Record the ice- edge is in the middle of both errors. Kriging employs a moving window 
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averaging procedure, which results in smoothing of steep edges. Consequently, the exact location of the ice 
edge cannot be inferred from altimetry, but the cliffs may be identifiable as a sequence of dense coverage. 
Comparing our map with the images in Wendler and others (1996), however, we note that the results of 
kriging altimeter data are surprisingly good. On their map in fig. 2, Mertz Glacier Tongue extends about 
80 km off the (averaged) coastline, Ninnis Glacier Tongue 20-30 km. (Also, 1993 - 85 = 8 years, 0.9 km per 
year for 1962-1993, 8x0.9 km =: 7.2 km or 8x1.2 km = 9.6 km advance since 1985/86, according to Wendler 
and others (1996) assuming constant rate of advance.) 

This indicates that, while SAR images are superior to altimetry-based maps for location of the ice edge, 
changes in advance/retreat can still be monitored from satellite-altimeter-derived digital terrain models, 
which are available for a longei time span (since 1978, Seasat mission). SAR imagery shows surface features 
in the floating ice tongue. The slope and elevation of the glacier surface are lost in the grey-shaded image. 
Additional information on ice lux can be derived from the surface elevation in the altimetry-based DTMs. 
calculated for the drainage basins. 

Another advantage of radar altimeter data over SAR data is that records have been available since the 1978 
Seasat mission. A maximal amount of information may be expected from a quantitative combination of SAR 
and radar altimeter data. 

In general, radar altimeter clcta and SAR data can be combined to study changes in position, surface 
morphology and elevation of Antarctic ice streams and glaciers. The advantage of RA data over SAR data 
is their longer record in time (since 1978), more frequent repeat, and lesser data volume, which facilitates 
frequent collection and rapid evaluation for large areas. 
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